
function [Hz,Hw,H,xs,DET1]=BlackBox2(Piecewise_Data,NH,ND,NP,v,Y,M,C,K)

X=Y(1:end-1);                                                              		% ���ۺϱ�������ȡ����Ӧ��Ϣ
w=Y(end);                                                                  		% ���ۺϱ�������ȡ��ת����Ϣ
%% ------------------A,Aw field start------------------------------
% construction of DET1,DET1w
     DET1=zeros(2*NH+1,2*NH+1);
     DET1w=zeros(2*NH+1,2*NH+1);

for ii=1:NH
   Temp_Matrix1=zeros(NH,NH);
   Temp_Matrix1(ii,ii)=1;
   det1=[0    -ii*w/v;
        ii*w/v    0;];
   det1w=[0   -ii/v;
          ii/v   0;];
   DET1(2:end,2:end)=DET1(2:end,2:end)+kron(Temp_Matrix1,det1);  
   DET1w(2:end,2:end)=DET1w(2:end,2:end)+kron(Temp_Matrix1,det1w);     
end

% construction of DET2,DET2w
    DET2=zeros(2*NH+1,2*NH+1);
    DET2w=zeros(2*NH+1,2*NH+1);

for ii=1:NH
   Temp_Matrix1=zeros(NH,NH);
   Temp_Matrix1(ii,ii)=1;
   det2=[ -(ii*w/v)^2   0;
        0          -(ii*w/v)^2;];
   det2w=[ -2*(ii/v)^2*w   0;
        0          -2*(ii/v)^2*w;];
   DET2(2:end,2:end)=DET2(2:end,2:end)+kron(Temp_Matrix1,det2);  
   DET2w(2:end,2:end)=DET2w(2:end,2:end)+kron(Temp_Matrix1,det2w);     
end

A=kron(DET2,M)+kron(DET1,C)+kron(eye(2*NH+1),K);
Aw=kron(DET2w,M)+kron(DET1w,C);
% ------------------A,Aw field end------------------------------
%% ------------------L field start------------------------------
Tau=2*pi*v;                                                                		% T=2*pi*v/w;
DTau=Tau/NP;                                                               		% DT=T/NP;
taus=(0:DTau:Tau-DTau)';                                                   		% ts=(0:DT:T-DT)';
L=zeros(NP*ND,(2*NH+1)*ND);
DL=zeros(NP*ND,(2*NH+1)*ND);
L(1:NP*ND,1:ND)=kron(eye(ND),ones(NP,1).*1/sqrt(2));
DL(1:NP*ND,1:ND)=kron(eye(ND),zeros(NP,1));
for jj=1:NH
   L(1:NP*ND,(2*jj-1)*ND+1:2*jj*ND)=kron(eye(ND),sin(jj*taus));            		% L(1:NP*ND,(2*jj-1)*ND+1:2*jj*ND)=kron(eye(ND),sin(jj*w/v*ts));
   L(1:NP*ND,2*jj*ND+1:(2*jj+1)*ND)=kron(eye(ND),cos(jj*taus));            		% L(1:NP*ND,2*jj*ND+1:(2*jj+1)*ND)=kron(eye(ND),cos(jj*w/v*ts));
end
% ------------------L field end------------------------------
    xs=L*X;                                                                		% Displacement
%% ++++++++++++++++++++++++ User Definition ++++++++++++++++++++++++++++++++
    gama=Piecewise_Data(7);                                                		% The setting of linear force
    fl = gama*cos(taus);		                                            	% linear force part
	
    k1=Piecewise_Data(1);k2=Piecewise_Data(2);B1=Piecewise_Data(5);B2=Piecewise_Data(6);
    fnl=zeros(length(taus)*ND,1);
    Dfnl=zeros(length(taus)*ND,1);	
	index1 = (xs>=-B1)&(xs<B1);													% Case 1
	index2 = (xs>=B1)&(xs<B2); 													% Case 2
	index3 = (xs>=-B2)&(xs<-B1);												% Case 3
	index4 = (xs<-B2);															% Case 4
	index5 = (xs>=B2);															% Case 5
	fnl = fnl - (k1.*double(index1).*xs );
	Dfnl = Dfnl - k1.*double(index1);
	fnl = fnl - (k1*B1.*double(index2) + k2.*double(index2).*(xs-B1));
	Dfnl = Dfnl - k2.*double(index2);	
	fnl = fnl - (-k1*B1.*double(index3) + k2.*double(index3).*(xs+B1));
	Dfnl = Dfnl - k2.*double(index3);	
	fnl = fnl - (k1*B1.*double(index4) + k2.*double(index4).*(B2-B1) + k1.*double(index4).*(xs-B2));
	Dfnl = Dfnl - k1.*double(index4);	
	fnl = fnl - (-k1*B1.*double(index5) + k2.*double(index5).*(-B2+B1) + k1.*double(index5).*(xs+B2));
	Dfnl = Dfnl - k1.*double(index5);	
    fs=fnl+fl;                                                             		% total forces
    FS=diag(Dfnl);     
    FV=zeros(length(xs),length(xs));

%  ++++++++++++++++++++++++    User Definition ++++++++++++++++++++++++++++++++++++++    
    BB=pinv(L)*fs;
    Hz=A-pinv(L)*FS*L-pinv(L)*FV*DL;
    Hw=Aw*X;
    H=A*X-BB;
end
